XID+

Posterior Predictive checking

In order to check whether our model seems ok, we can compare samples from the posterior joint distribution and use them to create 'different versions of the map' and compare them to the real map.

In [2]: 
field='COSMOS'
import pickle
import numpy as np
import pylab as plt
%matplotlib inline
obj=pickle.load(open('../XID_plus_output/XID+_'+field+'_obj.pkl', 'rb'))
In [3]: 
from scipy.sparse import coo_matrix
iter,chains,nparam=obj['chains'].shape
snsrc=obj['snsrc']
fit_d=obj['chains'].reshape(chains*iter,nparam)

Using JSanimation tool from Jake Vanderplas, we can create a movie of possible maps.

In [4]: 
x_range=np.max(obj['x_pix'])-np.min(obj['x_pix'])
y_range=np.max(obj['y_pix'])-np.min(obj['y_pix'])

im_temp=np.empty((x_range+1,y_range+1))
im_temp[obj['x_pix']-np.min(obj['x_pix']),obj['y_pix']-np.min(obj['y_pix'])]=obj['im_pix']
In [5]: 
from matplotlib import animation
from JSAnimation import IPython_display
In [6]: 
fig = plt.figure()
im = plt.imshow(im_temp,vmin=-5, vmax=100, cmap=plt.get_cmap('jet'))
plt.colorbar()
plt.xlim(20,200)
plt.ylim(20,200)
Out[6]:
(20, 200)
In [7]: 
def animate(i):    
    f=coo_matrix((fit_d[i,0:-1], (range(0,snsrc+1),np.zeros(snsrc+1))), shape=(snsrc+1, 1))
    rmap_temp=(obj['A']*f)
    x_range=np.max(obj['x_pix'])-np.min(obj['x_pix'])
    yrange=np.max(obj['y_pix'])-np.min(obj['y_pix'])

    im_temp=np.empty((x_range+1,yrange+1))
    im_temp[obj['x_pix']-np.min(obj['x_pix']),obj['y_pix']-np.min(obj['y_pix'])]=np.asarray(rmap_temp.todense()).reshape(-1)+np.random.randn(obj['snpix'])*obj['sig_pix']
    im.set_data(im_temp)
    return im,
In [8]: 
animation.FuncAnimation(fig, animate,
                        frames=200, interval=30)
Out[8]:


Once Loop Reflect
In [10]: 
plt.figure()
plt.hist(fit_d[:,-2])
Out[10]:
(array([  17.,   67.,  247.,  476.,  551.,  401.,  172.,   58.,    9.,    2.]),
 array([-2.84694539, -2.83449254, -2.82203969, -2.80958684, -2.79713399,
        -2.78468114, -2.77222829, -2.75977544, -2.74732259, -2.73486974,
        -2.72241689]),
 <a list of 10 Patch objects>)
In [19]: 
plt.scatter(fit_d[:,2761],fit_d[:,2762])
Out[19]:
<matplotlib.collections.PathCollection at 0x11361c8d0>
In [20]: 
hist,xedges,yedges = np.histogram2d(fit_d[:,2761],fit_d[:,2762],bins=40,range=[[0,50],[0,50]])
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1] ]
plt.imshow(hist.T,extent=extent,interpolation='nearest',origin='lower')
plt.colorbar()
Out[20]:
<matplotlib.colorbar.Colorbar instance at 0x113835200>
In []: 
print sx[]